#ifndef _WAVEFUN_COEFF_
#define _WAVEFUN_COEFF_


#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/filesystem.hpp>
#include <boost/mpi.hpp>
#include <boost/serialization/vector.hpp>
#include <cmath>
#include <cstddef>
#include <eigen3/Eigen/Dense>
#include <fstream>
#include <iostream>
#include <omp.h>
#include <stdlib.h>
#include <string>
#include <vector>
#include "../DF_Common/const_nucl.h"
#include "gsl/gsl_sf_gamma.h"
#include "gsl/gsl_sf_laguerre.h"
#include "gsl/gsl_sf_legendre.h"
#include "../DF_Common/progress_bar.h"
#include "../DF_Common/time_count.h"
#include "../DF_Mpi_Sep/info_wave.h"

//#include <boost/serialization/vector.hpp>
//#include "serial_control_wave.h"
//#include "serial_build_wave.h"

namespace mpi = boost::mpi;
namespace bf = boost::filesystem;

using Eigen::MatrixXd;
using std::string;
using std::vector;

struct nl_struct {
  int n12;
  int n3;
  int l12;
  int l3;
};

struct kpk_struct {
  int n_k;
  int n_kp;
};

struct k123_struct {
  int n_k12;
  int n_k3;
};

struct lbar_struct {
  int lbar;
  int lbar12;
  int lbar3;
};
struct klamb_struct {
  int lamb12_p;
  int lamb12;
  int lamb3_p;
  int lamb3;
};
struct S_lbar_struct {
  // [S][lbar] -> [slbar] for mpi
  int S_index;
  int lbar_index;
};
struct nlpnl_struct {
  // [nlp_num][nl_num] -> [nlpnl] for mpi
  int nl_p;
  int nl;
};

class Wave_Fun {

private:
  friend class boost::serialization::access;

  template <class Archive>
  void serialize(Archive &ar, const unsigned int version) {
    // ar &glbarS2_vec;
    // ar &glbarS2_mpi;
    // ar &mpi_layer;
    ar &mpi_size;
    ar &mpi_layer_PFK_vec;
    ar &PFKglbarS_vec_rank;
    ar &nlpnl2lay_rank;
  }

public:
  Wave_Fun(const int cut_lamb_t, const int hbar_omg_t, const int reg_pow_t,
           const int k_mesh_t, const int k_scale_t, const int p_mesh_t,
           const int Emax_final_t, const int Emax_init_t)
      : Cut_Lambda(cut_lamb_t), hbar_omega(hbar_omg_t),
        regulator_pow(reg_pow_t), k_mesh(k_mesh_t), k_scale(k_scale_t),
        p_mesh(p_mesh_t), Emax_final(Emax_final_t), Emax_init(Emax_init_t) {}

  const int k_mesh;
  const int k_scale;
  const int p_mesh;
  const int Emax_final;
  const int Emax_init;

  const double hbar_omega;    // = 22;
  const double regulator_pow; // = 2;
  const double Cut_Lambda;
  const double Mass_pi2 = Mass_pi * Mass_pi;

  Info_Wave IW;

  int lbar_max;
  int lbar_PFglbarS_max;
  int lbar_glbarS_max;
  int size_k123;
  int init_flag = -1;
  int evaluate_flag = -1;
  vector<vector<int>> E_Jac_conf;

  // mpi::communicator GlbarSCom;
  mpi::communicator WaveCom;
  int mpi_layer; // Don't changed initial in glbarS_vec_build_MPI and also used
                 // in PFKglbarS_build_MPI
  int mpi_size;
  vector<int> mpi_layer_PFK_vec; // was calculate PFKglbarS_build_MPI_II
  vector<vector<vector<double>>> glbarS2_mpi; // for mpi [layer][k123_p][k123]

  const vector<int> s1_vec{0, 2, 0, 2};
  const vector<int> s2_vec{0, 0, 2,
                           2};       // init at init() s1:[0,2,0,2] s2:[0,0,2,2]
  vector<S_lbar_struct> Slbar_state; // calculate in glbarS_init
  vector<vector<int>> Slbar_vec;
  vector<nlpnl_struct> nlpnl_state; // calculate in PFKglbarS_init
  vector<vector<int>> nlpnl_vec;

  vector<double> Pwave_w;                   // weight of intergration
  vector<vector<double>> Pwave_x;           // Pwave_x(0:2,0:k_mesh-1) x x^2 x^3
  vector<vector<vector<double>>> Pwave_vec; // P(n,l,1:k_max)
  vector<vector<vector<vector<double>>>> Pwave_vecwx; // P(0:2,n,l,1:k_mesh)

  vector<vector<double>> F_reg_vec;

  vector<vector<vector<vector<int>>>> nl_vec; // nl_vec[n12][n3][l12][l3] = num
  vector<nl_struct> nl_state;                 // nl_state[num].n12  .n3  .l3
  // glbarS_vec [S_index]   [lbar_state]      [kpk_state]    [kpk_state]
  //            s1, s2  lbar,lbar12,lbar3   k'_{12}k_{12}    k'_3 k_3

  vector<double> legendreP_x,
      legendreP_w; // will change in glbar_init and glbarS_init
  vector<vector<double>> legendre_Pl, legendre_PlW;
  // glbar_vec(0:1,0:lmax,0:k_mesh-1,0:k_mesh-1)
  vector<vector<vector<vector<double>>>> glbar_vec;
  // PFKglbar2_vec[num'][num][lambda1]['lbar']
  // lambda1=0,lbar=l3p ; lambda1=2, lbar=l3 :
  // PFKglbar2_vec[num'][num][lambda1][0] lambda1=1, lbar ={ max(l3p-1,l3-1) ,
  // min(l3p+1,l3+1)} : PFKglbar2_vec[num'][num][lambda1][lbar]
  vector<vector<vector<vector<double>>>> PFKglbar2_vec;
  vector<vector<double>> PFKglbar0_vec; //  PFKglbar0_vec[num'][num]

  // PFKglbarS_vec[num'][num][klamb_state][S_index][l_state]
  // --- IN 2pi - exchange term ---//
  vector<kpk_struct> kpk_state;   // kpk => {k'_{12}, k_{12}} or {k'_3, k_3}
  vector<k123_struct> k123_state; // k123 => {k_{12}, k_3 or {k'_{12}, k'_3}
  vector<vector<int>> k123_vec;   // k123_vec[k_12][k_3] = k123_state;
  vector<vector<double>> ksq;  // ksq[kpk_state][u] : (k'^2 + k^2 - 2k*k'*u)^0.5
  vector<vector<double>> ksq2; // ksq[kpk_state][u] : (k'^2 + k^2 - 2k*k'*u)^1
  vector<lbar_struct> lbar_state_g;       // lbar lbar12 lbar3
  vector<vector<vector<int>>> lbar_vec_g; // lbar_vec[lbar][lbar12][lbar3]
  // here is smaller that lbar_state & glbar_vec
  vector<lbar_struct> lbar_state;       // lbar lbar12 lbar3
  vector<vector<vector<int>>> lbar_vec; // lbar_vec[lbar][lbar12][lbar3]
  // [lbar][state]the power of k'_{12} k_{12} k'_3 k_3
  vector<vector<klamb_struct>> klamb_state;
  // klamb_vec[lbar][A][B][C][D]
  vector<vector<vector<vector<vector<int>>>>> klamb_vec;

  // PFKglbarS_vec[S_index][num'][num][lbar][klamb_state]
  // PFKglbarS_vec[S_index][nlpnl_index][lbar][klamb_state]
  vector<vector<vector<vector<double>>>> PFKglbarS_vec;
  vector<vector<vector<vector<vector<double>>>>>
      PFKglbarS_vec_T; // for test mpi
  // PFKglbarS_vec[mpi_layer][S_index][lbar][klamb_state]
  vector<vector<vector<vector<double>>>> PFKglbarS_vec_rank;
  vector<int> nlpnl2lay_rank; // calculated in PFKglbarS_build_MPI_II

  // glbarS_vec[S][lbar_state][kpk_12][kpk_3]
  vector<vector<vector<vector<double>>>> glbarS_vec;
  vector<vector<vector<vector<double>>>> glbarS_vec_T; // for test mpi
  // glbarS2_vec[S][lbar_state][k123_p][k123]

  vector<vector<vector<vector<double>>>> glbarS2_vec;
  vector<vector<vector<vector<double>>>> glbarS2_vec_T; // for test mpi

  //
  vector<vector<int>> Llbar_Ic_Reduced;
  vector<vector<int>> Llbar_Ic_Red_vec;

  double ho_k(int n, int l, double b, double k);
  /// HO. orbits in r space, radial part: R_nl(r)
  double ho_r(int n, int l, double b, double r);
  /// gauss-legendre quadrature
  void gauleg(const double x1, const double x2, vector<double> &x,
              vector<double> &w);
  //
  void init(int n_max, int l_max, int lbarS, int lbarS_G, int lbarS_P,
            int check_mpi, int MPI_TYPE);
  void init(int n_max, int l_max, int lbarS, int lbarS_G, int lbarS_P,
            int check_mpi, int MPI_TYPE, Info_Wave &IW_t);
  //
  void regulator_cal();
  //
  double glbar_cal(double k3, double k3p, int lbar, int g_index);

  //
  void PFKglbar0_build();
  void PFKglbar2_build();

  // void PFKglbarS_init_MPI5();
  //
  double PFKglbar0_cal2(int &n12p, int &n12, int &n3p, int &n3, int &l3p,
                        int &l3);

  double PFKglbar2_cal2(int &n12p, int &n12, int &n3p, int &n3, int &l3,
                        int &l3p, int lambda1, int &lbar);
  //
  // in 2-pi exchange term
  void glbarS_init(int lbar, int lbar12, int lbar3, int lmax_legendrePlW);
  double glbarS_cal(int s1, int s2, int lbar, int lbar12, int lbar3, int size_u,
                    int kpk_12, int kpk_3);

  void glbarS_vec_cal(int &S_index, int &lbar, int &lbar12, int &lbar3,
                      vector<vector<double>> &glbarS2_v);
  // calculate glbarS2_v[S_index][k123p_index][k123_index]
  void glbarS_vec_build();
  // build glbarS2_vec[S_index][l][k123p_index][k123_index]
  void glbarS_vec_build_MPI();    // Gather to one node
  void glbarS_vec_build_MPI_II(); // Gather to All nodes

  void PFKglbarS_init(int lbar, int lbar12, int lbar3, int Emax_F, int Emax_I);
  void PFKglbarS_init(int lbar, int lbar12, int lbar3, int Emax_F, int Emax_I,
                      vector<vector<int>> E_JC);

  double PFKglbarS_vec_cal2(int num_p, int num, int A, int B, int C, int D,
                            int S_index, int lbar, int lbar12, int lbar3);

  double PFKglbarS_vec_cal_II(vector<vector<double>> glbarS2_v, int num_p,
                              int num, int A, int B, int C, int D, int S_index,
                              int lbar, int lbar12, int lbar3);
  // glbarS2_v : glbarS_vec_cal(S_index, lbar, lbar12, lbar3, glbarS2_v);

  double PFKglbarS_vec_cal2_mpi(int num_p, int num, int A, int B, int C, int D,
                                int layer_num);
  double PFKglbarS_vec_cal_quick(vector<vector<double>> glbarS2_v, int &num_p,
                                 int &num, int &A, int &B, int &C, int &D,
                                 int &S_index, int &lbar_index);
  double PFKglbarS_vec_cal_quick(int flag, vector<vector<double>> glbarS2_v,
                                 int &num_p, int &num, int &A, int &B, int &C,
                                 int &D, int &S_index, int &lbar, int &lbar12,
                                 int &lbar3);
  void PFKglbarS_build();
  void PFKglbarS_build_MPI();
  void PFKglbarS_build_MPI_II();
  double PFKglbarS_quick(int &num_p, int &num, int &A, int &B, int &C, int &D,
                         int &S_index, int &lbar, int &lbar12, int &lbar3);

  double PFKglbarS_slow(int &num_p, int &num, int &A, int &B, int &C, int &D,
                        int &S_index, int &lbar, int &lbar12, int &lbar3) {
    vector<vector<double>> glbarS2_v;
    this->glbarS_vec_cal(S_index, lbar, lbar12, lbar3, glbarS2_v);
    int lbar_index = this->lbar_vec[lbar][lbar12][lbar3];
    return PFKglbarS_vec_cal_quick(glbarS2_v, num_p, num, A, B, C, D, S_index,
                                   lbar_index);
  };
  double PFKglbarS_rank(int &num_p, int &num, int &A, int &B, int &C, int &D,
                        int &S_index, int &lbar, int &lbar12, int &lbar3);

  double PFKglbarS_rank_quick(int &num_p, int &num, int &A, int &B, int &C,
                              int &D, int &S_index, int &lbar, int &lbar12,
                              int &lbar3);

  void print();
  void print(std::string filename);
  void print_Pmesh(std::string filename);
  /*
    void glbarS_save(string filename) {
      std::ofstream ofs(filename);
      // save data to archive
      boost::archive::binary_oarchive oa(ofs);
      // write class instance to archive
      oa << glbarS2_vec;
    };
    void glbarS_load(string filename) {
      std::ifstream ifs(filename);
      // save data to archive
      boost::archive::binary_iarchive ia(ifs);
      // write class instance to archive
      ia >> glbarS2_vec;
    };

    // void glbarS_save_file();
    // void glbarS_load_file();
    */
  void PFKglbarS_rank_save(string filename) {
    std::ofstream ofs(filename);
    // save data to archive
    boost::archive::binary_oarchive oa(ofs);
    // write class instance to archive
    oa << mpi_size;
    oa << mpi_layer_PFK_vec;
    oa << PFKglbarS_vec_rank;
    oa << nlpnl2lay_rank;
  };
  void PFKglbarS_rank_load(string filename) {
    std::ifstream ifs(filename);
    // save data to archive
    boost::archive::binary_iarchive ia(ifs);
    // write class instance to archive
    ia >> mpi_size;
    ia >> mpi_layer_PFK_vec;
    ia >> PFKglbarS_vec_rank;
    ia >> nlpnl2lay_rank;
  };

  void PFKglbarS_save_rank(int lamb_cut, int hbar_omg, int reg_pow, int k,
                           int ks, int p, int Emax_f, int Emax_i, int n, int l,
                           int lbarS, int lbarS_G, int lbarS_P);
  void PFKglbarS_save_rank(int lamb_cut, int hbar_omg, int reg_pow, int k,
                           int ks, int p, int Emax_f, int Emax_i, int n, int l,
                           int lbarS, int lbarS_G, int lbarS_P,
                           vector<vector<int>> E_JC);
  void PFKglbarS_load_rank(int lamb_cut, int hbar_omg, int reg_pow, int k,
                           int ks, int p, int Emax_f, int Emax_i, int n, int l,
                           int lbarS, int lbarS_G, int lbarS_P, int rank_t,
                           int mpi_size_t);
  void PFKglbarS_load_rank(int lamb_cut, int hbar_omg, int reg_pow, int k,
                           int ks, int p, int Emax_f, int Emax_i, int n, int l,
                           int lbarS, int lbarS_G, int lbarS_P, int rank_t,
                           int mpi_size_t, vector<vector<int>> E_JC);

  void glbarS_MPI_test();
  void PFKglbarS_MPI_test();

  //
  void evaluate();

private:
  int phase(int n); ///< return -1 if n is odd, 1 if n is even.
  int delta(int a, int b);
  double poww(double &a, int b); // quick power

  // matrix
  MatrixXd F_reg_mat;

  void glbar_init(int lmax);

  void legendre_compute_glr(int n, double x[], double w[]);
  void legendre_compute_glr0(int n, double *p, double *pp);
  void legendre_compute_glr1(int n, double x[], double w[]);
  void legendre_compute_glr2(double pn0, int n, double *x1, double *d1);
  void rescale(double a, double b, int n, double x[], double w[]);
  double rk2_leg(double t1, double t2, double x, int n);
  double ts_mult(double *u, double h, int n);
};

#endif
